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Bounded potentials are good models to represent the effective two-body interaction in some col- 
loidal systems, such as dilute solutions of polymer chains in good solvents. The simplest bounded 
potential is that of penetrable spheres, which takes a positive finite value if the two spheres are 
overlapped, being otherwise. Even in the one-dimensional case, the penetrable-rod model is far 
from trivial, since interactions are not restricted to nearest neighbors and so its exact solution is 
not known. In this paper the structural properties of one-dimensional penetrable rods are studied. 
We first derive the exact correlation functions of penetrable-rod fluids to second order in density at 
any temperature, as well as in the high-temperature and zero-temperature limits at any density. It 
is seen that, in contrast to what is generally believed, the Percus-Yevick equation does not yield 
the exact cavity function in the hard-rod limit. Next, two simple analytic theories are constructed: 
a high-temperature approximation based on the exact asymptotic behavior in the limit T — > oo and 
a low-temperature approximation inspired by the exact result in the opposite limit T — > 0. Finally, 
we perform Monte Carlo simulations for a wide range of temperatures and densities to assess the 
validity of both theories. It is found that they complement each other quite well, exhibiting a good 
agreement with the simulation data within their respective domains of applicability and becoming 
practically equivalent on the borderline of those domains. A comparison with numerical solutions 
of the Percus-Yevick and the hypernetted-chain approximations is also carried out. Finally, a per- 
spective on the extension of our two heuristic theories to the more realistic three-dimensional case 
is provided. 



I. INTRODUCTION 

In the last few years, the study of the structural and 
thermodynamic equilibrium properties of fluids with par- 
ticles interacting via strongly repulsive potentials have 
experienced a noticeable revival because of their interest 
in the physics of some colloidal systems. For instance, 
the effective interaction between two sterically stabilized 
colloidal particles can be accurately modeled by the hard- 
sphere potential^ 

On the other hand, the effective two-body interaction 
in other colloidal systems is much softer. For instance, 
the interaction potential for star polymers in good sol- 
vents can be shown to be ultrasoft, diverging only log- 
arithmically for short distancesi^ Furthermore, in the 
case of dilute solutions of polymer chains in good sol- 
vents, the centers of mass of two polymer chains can be 
separated by a distance smaller than the sum of their re- 
spective radii of gyration (and even can coincide at the 
same point), without violation of the excluded- volume 
conditions^ In such a case, the effective two-body po- 
tential is a bounded one, being well represented by the 
Gaussian core modeli^^iSiL& 

The simplest bounded potential is that of so-called 
penetrable spheres (PS), which is defined as 



where e > 0. This interaction potential was suggested 
by Marquest and Witten 9 as a simple theoretical ap- 
proach to the explanation of the experimentally observed 
crystallization of copolymer mesophases. In the last 
few years, the PS model has been the subject of sev- 
eral studies, 7 ! 10 ! 11 ! 12 ! 13 ! 14 ! 15 ! 16 ! 17 ! 18 ! 19 Density-functional 
theorjiiiii^ predicts a freezing transition to fee solid 
phases with multiply occupied lattice sites. The exis- 
tence of clusters of overlapped particles (or "clumps" ) in 
the PS crystal and glass was already pointed out by Klein 
et al.jiS who also performed Monte Carlo (MC) simula- 
tions on the system. In the fluid phase, the standard 
integral equation theories in general are not very reliable 
in describing the structure of the PS fluid, especially in- 
side the coreiii While the number of overlapped pairs is 
overestimated by the hypernetted-chain (HNC) theory, 
it is strongly underestimated by the Percus-Yevick (PY) 
theory. Other more sophisticated closures jiiii as well 
as Rosenfeld's fundamental-measure theory^ are able 
to predict the correlations functions with a much higher 
precision. Only in the combined high-temperature, high- 
density limit is the PS model amenable to an exact ana- 
lytical treatment^ 8 As applications of the model, let us 
mention that a mixture of colloids and non-interacting 
polymer coils, where the colloid-colloid interaction is as- 
sumed to be that of hard spheres and the colloid-polymer 
interaction is described by the PS model, has been 
studied^ The inhomogeneous structure of penetrable 
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spheres in a spherical pore has also been investigated^ 
In addition, some nonequilibrium properties have been 
analyzed recently^ and in the pastil 

As mentioned above, the classical integral equation 
theories (PY and HNC) do not describe satisfactorily 
well the structure of the PS fluid, especially inside the 
overlapping region, for the whole range of fluid densi- 
ties and temperatures. Thus, the PS model provides a 
stringent benchmark to test alternative theoriesii&iiii 
Even in the one-dimensional (ID) case, the PS model is 
far from trivial, since interactions are not restricted to 
nearest neighbors and so its exact solution is not known. 
A surprising consequence of the boundedness of the PS 
potential in the ID case is the plausible existence of a 
fluid-crystal phase transition, 18 thus providing one of the 
rare examples of phase transitions in ID systems^ 

Statistical mechanics has a long tradition of studying 
ID systems, especially in those cases where an exact so- 
lution to the many-body problem has been founds Of 
course, the ID PS model does not intend to describe all 
the properties of real polymers in solution, for which spa- 
tial dimensionality is known to be important.? 3 However, 
it seems worthwhile studying the ID PS model in order 
to understand some of the subtleties of the PS interac- 
tion and also to serve as a playground to test theoretical 
approaches that can be extended to the more realistic 3D 
case. 

The main aim of this work is to explore the possibility 
of constructing simple analytic theories for the structural 
properties of the ID PS fluid, based on known behaviors 
in the extreme cases of high and zero temperatures. Be- 
fore proposing those theories and in order to gain some 
insight, the exact properties in the limits of low density, 
high temperature, and zero temperature are worked out 
in Sec. |HJ some technical details being relegated to an 
Appendix. In particular, it is seen that the PY equa- 
tion does not yield the exact hard-rod correlation func- 
tions for all distances, in contrast to what is generally 
believediSiS^ As a simple extension to finite tempera- 
tures of the mean-field solution (asymptotically exact in 
the combined limit where the reduced temperature and 
density go to infinity, its ratio being kept finite) , we pro- 
pose a high-temperature (HT) approximation in Sec. lIIII 
A subtler task consists of the extension to finite temper- 
atures, imposing some basic continuity conditions, of the 
exact solution for hard rods. This is carried out in Sec. 
IIVI resulting in what we call the low-temperature (LT) 
approximation. Both approximations are compared with 
our own MC simulations in Sec. It is observed that 
both theories complement quite well each other, exhibit- 
ing a good agreement with the simulation data in their 
respective domains of applicability, which are wider than 
what one might have anticipated. To put the work on 
a broader perspective, we have also carried out compar- 
isons with numerical solutions of two integral equation 
theories, namely the PY and HNC theories. Both of them 
fail at low temperatures, but the HNC theory becomes 
very accurate and preferable to the PY theory at mod- 



erate and large temperatures. The paper ends in Sec. I VII 
with a discussion of the results and a perspective on the 
extension of the LT and HT approaches to the 3D case. 

II. EXACT PROPERTIES 

We consider in this paper a fluid of particles on a 
line interacting via the pairwise potential Hence- 
forth we take a = 1 as the length unit and define the 
reduced temperature as T* = ksT / 1. Even in this one- 
dimensional case, the exact solution of the problem for ar- 
bitrary density and temperature is not known, the main 
difficulty lying in the fact that one particle can inter- 
act simultaneously with an arbitrary number of parti- 
cles. Thus, it seems convenient to gain some insight by 
deriving a few exact results in limiting situations. 

A. Low-density limit 

Let us introduce the cavity function 

y(r) = e^/ k - T g(r), (2.1) 

where g{r) is the radial distribution function. In general, 
g(r) and y(r) depend parametrically on density (p) and 
temperature (T*). The virial expansion of the cavity 
function reads 

oo 

y(r) = l + Y,yn(r)p n , (2.2) 

n=l 

where the coefficients y n {r) depend parametrically on T* 
and are represented by diagrams In particular, 




(2.4) 

Here, the open circles represent root points separated 
by a distance r, the filled circles represent field points 
to be integrated out, and each bond represents a Mayer 
function 

/(r) = e -p(r)/k B T _ L ( 2 5 ) 

Thus, for instance, 

lt /\ = / dr 3 /(r 13 )/(r 23 ), (2.6) 
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Label Diagram 



< r < 1 



1 < r < 2 



2 < r < 3 




= / dr 3 / dr 4 /(ri3)/(r34)/(r 2 4)/(ri4). 



(2.7) 

So far, Eqs. (|2.1() - (|2.7JI hold for any interaction po- 
tential and any dimensionality. In the special case of 
penetrable spheres, Eq. (|2.1|l yields 



l(r) 



(l-x)y(r), r < 1, 



K v(r), r>l, 
where we have called 

Moreover, the Mayer function becomes 
f(r) = xfus(r), 

where 



/ H s(r) 



-1, 



0. 



r < 1, 
r > 1, 



(2.8) 

(2.9) 
(2.10) 

(2.11) 



is the Mayer function of hard spheres. Therefore, the spa- 
tial dependence of each one of the diagrams contributing 
to the virial expansion (|2.2|) is exactly the same as for 
hard spheres. The only difference is that each diagram is 
now multiplied by the temperature-dependent parameter 
x raised to a power equal to the number of bonds in that 
particular diagram. 

In the special case of one-dimensional penetrable 
spheres (i.e., penetrable rods), it is not difficult to eval- 
uate the integrals represented by the diagrams in Eqs. 

and • The results are displayed in Table |T| 
Note that diagram 2C is the square of diagram 1. All 
the functions are continuous everywhere and vanish for 
r > 3. At r = 1, diagrams 2A, 2B, and 2D have a 
second-order discontinuity. At r = 2, diagrams 1 and 2B 
have a first-order discontinuity, while 2C and 2D have 
a second-order discontinuity. Finally, 2A has a second- 
order discontinuity at r = 3. 

Inserting the expressions displayed in Table[I]into Eqs. 
(12. 3|) and l|2.4|) we get the exact results shown in Table 
ITT1 The discontinuities at r = 1, 2, and 3 are 

y"(l + ) - y"(l-) = -p 2 x 3 (3 -4x + x 2 ) + 0(p 3 ), (2.12) 



y'(2 + )-y'(2-)=px 2 + 0(p 2 ), (2.13) 



y"(3+)-y"(3-)=A 3 + 0(p 3 ). (2.14) 

It is interesting to compare the exact density expansion 
with the HNC and PY approximationsi2£i2I In the HNC 




x 2 (2-r) x 2 (2-r) 



2A 



2B 



2C 



2D 



-x J 3 - r 2 



1(3 - r f -^(3-r) 2 




%-(6-2r-r 2 ) 2L- (2 - r) (4 - r) 



x 4 (2-r) 2 x i (2-r) 2 



-x 5 (3-2r) -x 5 (2-r)' 2 



TABLE I: Explicit expressions of the diagrams contributing 
to j/i (r) and y 2 (r) . 



approximation, the so-called bridge (or elementary) dia- 
grams are neglected. To second order in density, the only 
bridge diagram in Table[I]is 2D. Therefore, the function 
2/1 (r) is retained but the function j/aC?") is approximated 
by the function y™ (r) given in Table ILTI In the PY ap- 
proximation, apart from the bridge diagrams, a subset of 
the remaining diagrams is also neglected. In particular, 
the PY expression for j/2 (r) neglects diagrams 2C and 2D 
of Tabled what results in the function y^ Y (r) shown in 
Table CH 

In the hard-rod limit (e — > oo or, equivalently, T* — > 
or x — > 1), the exact, HNC, and PY functions y%{r) re- 
duce to the expressions also included in Table [H] While 
the PY approximation for hard rods yields the exact pair 
correlation function for r > 1 (and this happens to any 
order in density because of a mutual cancelation of the 
neglected diagrams), we see that the PY cavity function 
fails inside the overlapping region (r < 1). In principle, 
this is not a serious drawback for strict hard rods, since 
g(r) = for r < 1, regardless of the expression of the 
cavity function. However, if we consider penetrable rods 
at low temperature (T* <C 1 or 1 — x <C 1), then g(r) for 
r < 1 takes on relevant (although small) non-zero values, 
so that the expression of lim^.^o U2 Y ( r ) f° r f < 1 antic- 
ipates the sort of limitations one can expect from the so- 
lution of the PY theory to penetrable-sphere modelsiiii^ 
It is worth noting that the second-order discontinuities 
at r = 1 of diagrams 2A, 2B, and 2D of Table [Q exactly 
compensate each other when those diagrams are added 
to get limy* U2 (r). On the other hand, since diagram 
2D is neglected in the HNC and PY approximations, an 
artificial second-order discontinuity at r = 1 remains in 
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Function 



< r < 1 



1 < r < 2 



2 < r < 3 



2/1 M 
1/2 ( r ) 

y 2 HNC M 

lim T »_ y 2 (r) 
Um T »_ 1/J Y (r) 



a; 2 (2 - r) 

-x 3 (3-r 2 ) +x 4 (8-4r- ir 2 ) - ( { 
-x 3 (3 -r 2 ) + x 4 (8 -4r - ^r 2 ) ' 
-I 3 (3 - r 2 ) + x l (6 - 2r ~ r 2 ) 
ir 2 -3r+ J 
ir 2 - 4r + 5 
3 - 2r 



a; 2 (2-r) 

'(r - 3) 2 + x 4 (2 - r) (5 - |r) - # (r - 2) 2 -a£(r - 3) 2 

_4( r -3) 2 +* 4 (2-r)(5-§r) -#(r-3) 2 

-#(r- 3) 2 +x i (2 - r)(4 - r) _^( r - 3 ) 2 

Ir 2 -3r+I -^(r--3) 2 

r 2_ 5l . + u _i( r _ 3 )2 

|r 2 -3r+| -^(»--3) 2 



TABLE II: Explicit expressions of the functions yi(r) and y2{ r )- The first, second, and fifth rows are exact results, while the 
third, fourth, sixth, and seventh rows correspond to the HNC and PY approximations. 



those approximations. It is interesting to point out that, 
although the hard-rod interaction (e — > oo) is more sin- 
gular than the penetrable-rod interaction (e = finite), 
the cavity function y(r) of the former is paradoxically 
analytic at r — 1, while that of the latter presents a 
second-order discontinuity. When considering the radial 
distribution function g{r), however, a jump discontinuity 
exists in both cases. 



0.02 



5^ 



0.00 



-0.02 



-0.04 . 



T=2 



exact 
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r 

FIG. 1: (Color online) Plot of the virial coefficient j/2(r) at 
T* = 2 (top panel) and T* — (bottom panel). The solid 
lines are the exact results, the dashed lines are the HNC pre- 
dictions, and the dotted lines are the PY predictions. 

Figure ^ compares the exact function j/ 2 ( r ) with the 
HNC and PY approximations at T* = 2 and T* = 0. 
Although restricted to low densities and one-dimensional 



systems, Fig.ndearly illustrates some of the general fea- 
tures found in three-dimensional penetrable spheres at fi- 
nite densities ii&ii the PY underestimates the penetrabil- 
ity effect (r < 1), while the HNC approximation overesti- 
mates it. Moreover, the HNC theory tends to be better at 
higher temperatures, while the PY is preferable at lower 
temperatures. If we characterize the quality of each ap- 
proximation by the separation of the corresponding 7/2 (0) 
from the exact result, it turns out that the temperature 
beyond which the HNC approximation becomes better 
than the PY approximation is T* = l/ln3 ~ 0.91. If we 
take instead y 2 (l) as the characteristic quantity, then the 
preceding temperature is replaced by T* = 1/ In 2 ~ 1.44. 
In either case, T* « 1 seems to be the typical tempera- 
ture beyond which the HNC approximation prevails over 
the PY one. This will be confirmed later in Sec. |Vj 

Although in this paper we are more interested in the 
structural properties than in the thermodynamic ones, 
it is worthwhile obtaining the first four virial coefficients 
as functions of temperature. For penetrable rods, the 
compressibility factor Z, the isothermal susceptibility \i 
and the excess internal energy per particle related 
to the cavity function as follows: 



p 



l + pxy(l), 



(2.15) 



dp 
dp 



= 1 + 2 p[J q dr[(l - x)y{r) - 1] 



+ J dr[y{r) -1}}, 



u cx = pe(l -x) dry(r), 
Jo 



(2.16) 



(2.17) 



where p is the pressure and f3 = l/fc^T. The virial coef- 
ficients B n (x) arc defined by the density expansion 



Z = l + Y j B n {x)p n - x . 



(2.18) 



n=2 



Inserting the exact cavity function to second order in 
density into Eqs. H2.15JI - H2.17JI . an d taking into ac- 
count the thermodynamic relations x _1 = 9{pZ)/dp and 
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Theory 
Exact 
PY-v 
PY-c 
PY-e 
HNC-v, HNC-e 
HNC-c 



Bj(x) 



(2 - 3x) 



4 _ r 

3 3 
4 Ir. 14 



-a; 4 (2 



-a 4 (2 -I*) 



gj(l) 
1 
1 
1 

4 
■ r i 

3 
2 

U 

12 



To 
1.0120 
0.9102 
1.1802 
0.7982 
1.1802 
0.8640 



1.4427 
1.3121 
1.6369 
1.1802 
1.6369 
1.2573 



BA . 

I mm 

-0.02344 
-0.03236 
-0.01165 
-0.04265 
-0.01747 
-0.03622 



TABLE III: Fourth virial coefficient B±(x) and other related quantities as given exactly and by the PY and HNC approximations 
through the virial (v), compressibility (c), and energy (e) routes. 



pdu C x/dp — dZ/d{3, one consistently gets B 2 {x) 



x, 



B 3 (x) 



— „,3 



and the fourth virial coefficient B^{x) dis- 



played in Table 11111 On the other hand, insertion of 
y PY (r) yields a different B^x) depending on whether 
the virial (v) route IJ2.15JI . the compressibility (c) route 
(|2.16(l . or the energy (e) route 12.17|l is used. This 
internal inconsistency is also present in the case of 
y H (r), except that now the virial and energy routes 
are equivalent^ These approximate expressions for the 
fourth virial coefficient are also included in Tablc lTTTl The 
exact and the approximate B^(x) share the properties of 
changing sign at a certain temperature T * and having 
a (negative) minimum value -B 4 | min at a higher temper- 
ature Tj* in . Only the PY-v and PY-c approximations 
are consistent with the exact hard-rod value B^l) = 1, 
while only the PY-c approximation fails to yield the cor- 
rect limit B4 — > —2a; 4 for high temperatures. From Table 
IIIII one can conclude that the best global agreement with 
the exact B^(x) is achieved by the PY-v approximation. 
It is interesting to note that, while the energy route to the 
equation of state is ill defined for hard spheres, one can 
circumvent this problem by first getting B^{x) through 
the energy route and then taking the zero-temperature 
limit lim x ^!B 4 (a;). The fact that B% Y - e (l) ^ B FY - V (1) 
shows that the equivalence between the energy and virial 
routes proven in Ref. |2^ when hard spheres are obtained 
from the square-shoulder potential does not hold when 
hard spheres are obtained from penetrable spheres. 



B. High-temperature limit 

In the high-temperature limit (T* — > oo), the parame- 
ter defined by Eq. I|2.9|l tends to zero, i.e., x w 1/T* — > 0. 
Since a diagram having m bonds is of order x m , we can 
neglect all the diagrams contributing to a given coeffi- 
cient except the one having the least number of 
bonds (namely, m = n + 1 bonds). In other words, only 
the linear chain diagrams survives 1 ^ 



Vn{r 12) 



1 3 4 
o • • 



n+2 2 
— • O 



x n+1 J dr 3 J dr 4 ...y dr„ +2 

X/Hs(ri3)/Hs(^34) ' ' ■ / H s(?n+2,2)(2.19) 



In Fourier space, 



y n (k)^x n+L / HS (fc) 



n+l 



(2.20) 



where /hs(^) is the Fourier transform of the hard-sphere 
Mayer function (|2~H1) . From Eqs. lO and (|2~2U|) it is 

straightforward to get the Fourier transform of y(r) — 1, 
so that we finally have 

y{r) -> 1 + xw(r), (2.21) 

where w(r) is the inverse Fourier transform of 

[/hsW] 2 



w(k) = px 



1 - pxfm(k) 



(2.22) 



The limit result (|2.21l) can be written in a number 
of equivalent ways. For instance, the total correlation 
function h(r) = g(r) — 1 becomes, on account of Eqs. 
(|2~%|) . J5TTJ, and'(E2TJ), 



h(r) -» x [w(r) + fns(r)} 
In Fourier space, 



h(k) — > x 



1 - pxfm{k) 



(2.23) 



(2.24) 



where use has been made of Eq. (|2.22|) . From Eq. (|2.24|1 it 

is straightforward to get the structure factor S(k) = 1 + 
ph{k), as well as the Fourier transform c(k) = h(k) /S(k): 



S(k) 



1 - pxfm(k) 



c(k) -» z/hs(*0. (2-25) 



The last expression in Eq. (|2.25l) is equivalent to c(r) — > 
fir)- 

The asymptotic behaviors H2.21JI - H2.25J) are of mean- 
field typo 7 i 30 i 31 i 32 i 33 and hold in the combined limit T* — > 
00, p — > 00 with px ~ p/T* — fixed for any dimensional- 
ity. In the one-dimensional case, one simply has 



f m (k) = -2 sink /k, 



(2.26) 
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so that 
w(r) = — 

7T 



dk coskrw(k), w(k) = 



4px sin k 



- 2pxk sin k 
(2.27) 

It is proven in the Appendix that the density expansion 
of w(r) is 



where 

w n {r) - 



E 



(-1) 



n=2 



n-\-m 



(2.28) 



— ' m\(n — m)) 

m=0 v ' 



(?i-2m-r) n " 1 e(fi-2ra-r), 
(2.29) 

0(a:) being the Heaviside step function. It is also shown 
that 



>W>n(l) 



Jn+1 



(0) 



1 



n > 2, 



(2.30) 



this property being needed to prove the thermodynamic 
consistency between the virial and energy routes to the 
equation of stated 

The virial series (I2.28|l converges for px < v{r), where 
l/v(r) = linin^oo \w n +i{r) / w n (r)\. Equation (|2.30|) im- 
plies that i/(0) = v(l), but otherwise the radius of con- 
vergence v(r) could in principle be r-dependent. We 
have numerically checked that this is not the case and 
that v{r) — \ regardless of the value of r. From Eq. 
(|2.27|l we see that the radius of convergence v(r) = h 
is a consequence of the mathematical singularity of w(r) 
at the negative density 2px = — 1. Therefore, we con- 
clude that in the high-temperature domain the virial ex- 
pansion of the cavity function converges uniformly for 
p < l/2x w T* /2. 

Using Eqs. (|2~2"T|) and (|2~2~g|) in Eq. (|2~Tojl . we obtain 
the virial expansion of the equation of state in the high- 
temperature limit, 



P 



pksT 



1 



px - 



oo 

E 

n — 3 



p 



X W r , 



l(1), 



(2.31) 



so that the virial coefficients are Bi (x) = x and B n (x) — » 
x n w n -i(l) for n > 3. Therefore, the virial series of the 
equation of state converges for p < (2x) _1 in the high- 
temperature limit. 

By an adequate reordering of terms (see the Ap- 
pendix), the function w(r) can be written in real space 
as 



;(r) = 6(1 - r) + 2 ^(-l) p J' p (2px(p - r))0(p - r), 

P =i 

(2.32) 

where J' p {z) is the first derivative of the Bessel function 
of the first kind J p (z)^4 In particular, in the overlapping 
region (r < 1), we have 

oo 

«;(r < 1) = 1 + 2 J' P (2px(p - r)). (2.33) 

P =i 



Likewise, in the shell n < r < n + 1 with n > 1, 

oo 

w(n< r < n + l) = 2 ^ (-l) p J p (2px(p - r)). (2.34) 

p— n+l 

Using known properties of the Bessel functions it is pos- 
sible to derive some exact results. For instance, the 
second-order, first-order, and second-order discontinu- 
ities of w(r) at r = 1, 2, and 3, respectively, are 

to"(l+)-tfl"(l-) = -i{pxf, 



</{2 + )-w'{2-) = px, 



(2.35) 



(2.36) 



to" (3+) - to" (3~) = (/m:)' 



(2.37) 

Equations l|2~3l>|) ~ (j2~3?|) give the singularities of the cav- 
ity function for any density in the high-temperature limit. 
Of course, they are compatible with Eqs. (j2.12|l - (|2.14(l 
for low densities. Other interesting results are 

w'(1)=/iiK0)-ib(2)-1], (2.38) 
to' (2+) = px[w(l) - to(3)] , (2.39) 
w'{n) = px [w(n - 1) - w(n + 1)] , n > 3. (2.40) 
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*> v • 
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FIG. 2: (Color online) Plot of the function w(r) at px = 0.3, 
0.5, and 0.8. 

Although the series representation (|2.32|l is useful to 
derive explicit results and seems to converge for any value 
of px, it is more practical to use the integral representa- 
tion H2.27f) from a computational point of view. Figure 
121 shows the function to(r) for px = 0.3, px — 0.5, and 
px = 0.8. As expected, the function w(r) becomes more 
structured as the product px grows. 
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C. Zero-temperature limit: Hard rods 

The correlation functions for one-dimensional fluids 
with interactions restricted to nearest neighbors are ex- 
actly knowni 2 ^ 3 ^ 3 9i 3 7i 38 i 3 ? However, as indicated in Sec. 
[I] the bounded nature of the potential allows for 

multiple pair interactions and, as a consequence, the in- 
teractions are obviously not restricted to nearest neigh- 
bors. 

On the other hand, as the temperature is lowered, the 
occurrence of simultaneous overlapping of three or more 
particles becomes less and less likely and the nearest- 
neighbor approximation becomes more and more accu- 
rate, until the exact result is obtained at zero tempera- 
ture (hard-rod limit). In that limit (T* = or x = 1), 
the radial distribution function g(r) vanishes inside the 
core (r < 1) but the cavity function y(r) does not and 
so special care must be taken to get lim-.^o Z/( r ) as a 
function of density for r < 1. As already recognized by 
Stell^S in order to get the hard-sphere cavity function for 
r < 1 one must take the limit of a suitably parameterized 
non-singular potential. 

We start by recalling the exact solution in the case 
of nearest-neighbor interactions. At arbitrary density 
and temperature, the one-dimensional radial distribution 
function can be written as 38 ' 39 



9{r) 



1 X 



(2.41) 



n=l 



where p n (r) is the probability density distribution that 
the n-th neighbor of a given particle is a distance r apart. 
The distributions {p~(r)} obey the recursive relation 

Pn(r)= [ dr'p n -x{r') Pl (r-r'), n>2. (2.42) 



The explicit form of the nearest-neighbor probability dis- 
tribution is 



Pi 



(r) = Ke-^/ kBT e~ ir , 



(2.43) 



where the amplitude K and the damping constant £ will 
be determined later on. 

The convolution in l|2.42|) suggests the introduction of 
the Laplace transforms 

/>oo />oo 

P n {t)= dre- rt p n (r), G{t) = dr e - rt g(r). 
Jo Jo 
(2.44) 

Thus, Eq. (|2~4'3|) yields 

P n (t) = Pn-i^P^t) = P{ l (t), n>2. (2.45) 
When this is introduced into Eq. (|2.41|l one gets 



ow-i .ijM) 

n— 1 



(2.46) 



From Eq. I|2.43|l . it is obtained 

i\(f)=JHi(t + 0> 



where 



fi(t) = / dre~ rt e-^ r ^ kBT 
Jo 



(2.47) 



(2.48) 



is the Laplace transform of exp[— <£>(")/fc-T]. 

To close the solution, one needs to determine the pa- 
rameters K and £. This is done by imposing basic consis- 
tency conditions. The relationship between G(t) and the 
Laplace transform H(t) of the total correlation function 
h(r) = g{r) - 1 is 



G(t) = - t +H(t). 



(2.49) 



Since H(t) must be finite at t = (from the compress- 
ibility route to the equation of state), it follows that 

G(t) = \ + O{t ), P 1 {t) = l--t + 0(t 2 ). (2.50) 
t p 

The last equality in (I2.5U|) implies that Pi(0) = 1 and 
P[(Q) = -1/p. Hence, K = and £ is tne solution 

of the equation 



(2.51) 



We emphasize that Eqs. H2.41j l~ H2.51j l provide the ex- 
act solution only in the case of nearest-neighbor inter- 
actions. Let us now "forget" for the moment that the 
penetrable-sphere potential (f 1 . 1 1> is not restricted to near- 
est neighbors and apply the above scheme to it. In that 
approximation, Eqs. (j2.43j) and (j2.48j) yield 

(2.52) 



Pl (r) = Ke-^ r [l-x + xQ{r - 1)] , 



n(t) = 



1 — x + xe 



(2.53) 



The damping coefficient £ is the solution of the transcen- 
dental equation H2.51jl . i.e., 



(£-0 = 



1 -x 

where we have called 



£' = - - 1. 

P 

In addition, the amplitude K becomes 



K 



and Eq. (|2~47jl gives 
Px{t) 



i + £ 



(2.54) 



(2.55) 



(2.56) 



(2.57) 
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Inserting this into (|2.46(l and expanding in powers of e 
we get 



P G(t) = 



7 + 



1 



n+l 



Its Laplace inversion yields 

oo 



where 



n=0 



/« p -('r r n-l 



C e 



P 



[n + (£ - Or] 



(2.58) 



(2.59) 



(2.60) 



In particular, using Eq. I|2.8|) . the cavity function inside 
the core is, in this nearest-neighbor approximation, 



/v ' 1-x p(l-x) 



(2.61) 



Of course, Eqs. (|2.52|I - H2.61|I are not exact at finite 
temperature. However, they become exact in the hard- 
rod limit (T* — * or x — » 1). In that limit, the solution 
of Eq. f231jl is f = $o - Co e?0 (1 - a;) + • • • , so that £' = 
£ — Coe^ n (1 — i) + • • • , where 



Therefore, 

lim tp n (r) 



so 



/0 



1-P 



e -for r n-l 



(2.62) 



n > 1, (2.63) 



lim y(r < 1) = — e 



£o .-$ (r-l) 



(2.64) 



The final explicit expression of the cavity function for 
hard rods is then 



l im y M = » e -Co(r-l) + V- *0 -fo(r 
r*^o v p ^ p(n - 1)! 



xir-n^Qir-n). 



(2.65) 



This equation not only gives the exact hard-rod radial 
distribution function g{r) = y(r)Q(r — 1), but also the 
cavity function inside the hard core (r < 1). Interestingly 
enough, y(r) for r < 1 is just the analytical continuation 
of its expression for 1 < r < 2, so that y(r) is analytical 
at r = 1, even though the potential is highly singular 
at that point. This singularity manifests itself at r = n 
with n > 2, where y(r) = g(r) presents a discontinuity of 
order n — 1. 

We are not aware of any previous derivation of Eq. 
By expanding in powers of density, it is easy 



to verify that Eq. (|2.65|l is consistent with the exact 
liniT*— >o Ui( r ) an d hmy*^oy2( r ) given in Table ITU As 
a further test of the exact character of Eq. Ij2.64|) . let us 
check the fulfillment of the following two zero-separation 
theorems^ as applied to hard rods, 



In 2/(0) = 



dlny(r) 



Or 



(2.66) 



(2.67) 



In Eq. (|2.66() . /i cx is the excess chemical potential, which 
is given by 

-i-p cx = Z{p) - 1 + T &p' W ~ l , (2.68) 
kbT J p' 

where Z = l+py(l) is the hard-rod compressibility factor 
[cf E q. f2~T3|l ]. Since y(l) = 1/(1 - p), Eqs. (|2~^|l and 
(|2~B7|) imply that 



y(0) = - =»(l)eWW, 



dhiy(r) 



Or 



r=0 



1 



(2.69) 



(2.70) 



It is straightforward to check that Eq. I|2.64|l verifies both 
(l2~?)51) and JOB- 

In a general one-dimensional system, the Fourier trans- 
forms of the total and direct correlation functions can be 
obtained from the Laplace transform of the radial distri- 
bution function as 



h(k) = G(t = ik) + G(t = -ik), c{k) 



h{k) 



1 



ph(k) 
(2.71) 

Inserting Eq. I|2.57|l into Eq. I|2.46|l and taking the zero- 
temperature limit one finally gets 



lim c(k) = 

p 



1 — cos k sin k 



lim c(r) = 

T>->0 



k 2 



1 — pr 



9(1 -r). 



(2.72) 



(2.73) 



The bridge function is defined a.a£lB(r) = In y(r)+c(r) — 
h(r). Therefore, Eqs. l(2~Mjl and (|2~73l yield 

lim B(r < 1) - Q - ln(l - p). (2.74) 

T*^0 (1 - p) 2 

As is well known, the PY theory yields the exact corre- 
lation functions g(r) and c(r) for hard rodsi^ii2S42 How- 
ever, it provides a wrong cavity function y(r) inside the 
hard core, as already seen to order p 2 in Sec. Ill Al To 
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show this for finite densities, we recall that the PY ap- 
proximation consists of the closure^! y(r) = g(r) — c(r). 
Therefore, 



lim y FY (r < 1) 



-c(r < 1) 



1 



pr 



(I-P) 2 



(2.75) 



This result, which gives rise to a second-order disconti- 
nuity at r = 1, differs markedly from the exact result 
(|2.ti4|) . The PY cavity function lim^.^o U PY ( r ) verifies 
the zero-separation theorem H2.70II but not H2.69(l . As for 
the PY bridge function inside the core, the result is 



hm B (r < 1) = ^ — 

t*^o (1 - p) 2 

which differs from Eq. (|2.74(l . 

100 



X 

c 

T 



•ln 



pr 



, (2.76) 




Before closing this Section it is worthwhile noting a 
flaw in the "proof" given in Ref.|25|on the exact character 
of Eq. (|2.75|) . It is first proven that the PY closure is 
consistent with the exact hard-core radial distribution 
function or, equivalently, with the exact cavity function 
for r > 1. Next, taking the limit r — > 1 + and taking into 
account the continuity of y{r) at r = 1, it is argued that 



1 — pr 
Jl^pf 



= lim lim y(r) = lim 
r^l+T*^Q r->l- 

= lim y(r < 1). 



lim y(r) 



(2.77) 



While the two first equalities are entirely correct, the 
third one is not justified in general, so that the conclusion 
that lim T .^ y( r < 1) — lim T »^ y PY { r < 1) is flawed. 



III. HIGH-TEMPERATURE APPROXIMATION 

In the limit of asymptotically high temperatures, the 
cavity function becomes the mean-field result l|2.21|l . 
where the function w(r), which in the one-dimensional 
case is given by Eqs. (I2.27|) or l|2.32ll . depends on density 
and temperature through the product px only. It is then 
natural to expect that the simple ansatz 



y(r) — 1 + xw(r) 



(3.1) 



constitutes a good approximation for sufficiently large 
temperatures. As a matter of fact, comparison with 
three-dimensional simulation results shows that this in- 
deed the case for T* > 3i 

The approximation neglects all the diagrams dif- 
ferent from the linear chain ones. A better approximation 
can be expected by making 



y(r) = exp[xw(r)], 



(3.2) 



which reduces to l|2.21|l in the limit x — + 0. In addition 
to the linear chain diagrams, the approximation (|3.2|) re- 
tains all the reducible diagrams that factorize into linear 
chain diagrams. In particular, l|3.2|l implies that 



FIG. 3: (Color online) Plot of the cavity function y(r) and 
the bridge function B(r) in the hard-rod limit at p — 0.8. The 
solid lines are the exact results, while the dashed lines are the 
PY predictions. 

Figure shows limT*^o y(r) and liniT*^o i?(?™) at p = 
0.8. It clearly illustrates the tendency of the PY ap- 
proximation to underestimate the cavity function inside 
the core, a feature that continues to hold in the case of 
three-dimensional penetrable spheres^ It is also quite 
apparent that the HNC closure B(r) =0 is far from be- 
ing accurate in the hard-rod limit, especially for r < 1. 



2/2 (r) 




(3.3) 













O ( 


+ 

J c 






n 



2/3 (r) 

A third possibility is the Pade approximant 

1 



1 — xw(r) 



(3.4) 



(3.5) 
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In that case, the diagrams retained are the same as in 
the approximation l|3.2[) . except that the numerical fac- 
tors in front of the reducible diagrams are larger. For 
instance, the factors affecting the second diagram of Eq. 
(|3.3() and the second and third diagrams of Eq. (|3.4II are 
1 (instead of |), 2 (instead of 1), and 1 (instead of i), 
respectively. This can be interpreted as an effective way 
of compensating for the neglected diagrams. 



(|3.2() and the linear approximation l|3.1|) . Therefore, in 
what remains we will adopt Eq. (|3.5|) . supplemented by 
Eq. I|2.27|l . as our high-temperature (HT) approximation. 
The highest value of w(r) occurs at r = and is an in- 
creasing function of px. For px < 1.1 one has w(0) < 1 
and so Eq. (|3.5|) is mathematically well defined for any 
temperature if px < 1.1. This is the range where the 
fluid is expected to be stable against the solid^ 
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r 

FIG. 4: (Color online) Plot of the cavity function y(r) at p = 3 
and T* = 10 (top panel) and T* = 5 (bottom panel). The 
solid lines are results obtained from MC simulations, the dot- 
ted lines represent the approximation 13. H . the dashed lines 
represent the approximation 13.21 . and the dashed-dotted 
lines represent the approximation 13.51 . 

The three heuristic approximations l|3.1|l . I|3.2|l . and 
(|3.5|) . among other possible ones, have the potential of 
describing satisfactorily well the structure of penetrable- 
rod fluids at sufficiently high temperatures. To test this 
expectation and also to determine which approximation 
is the most accurate one, we have compared them with 
MC simulations at T* = 10 and T* = 5 in Fig. H It can 
be observed that the three approximations tend to under- 
estimate the correlation function for r < 2, this effect be- 
coming more pronounced as the temperature decreases. 
On the other hand, the Pade approximant l|3.5|l is closer 
to the MC curves than the exponential approximation 



IV. LOW-TEMPERATURE APPROXIMATION 
A. Construction of the approximation 

We have already mentioned that the automatic trans- 
lation to penetrable rods [cf. Eqs. 12.52fl - H2.61l) ] of the 
scheme valid for nearest-neighbor interactions does not 
yield the exact correlation functions, except at zero tem- 
perature. On the other hand, one might reasonably won- 
der whether Eqs. H2.52fl - I|2.61|) constitute a good approx- 
imation for (non-zero) low temperatures (i.e., T* C 1 or 
1 - i < 1). A drawback of the approximation (|2.52(l - 
(|2.rill) . however, is that it yields a cavity function that 
has a jump at r = 1, i.e., y(l + ) ^ y(l~), while the exact 
y{r) and its first derivative must be continuous at r = 1. 

The failure to satisfy the condition y(l + ) = 2/(l~) can 
be traced back to the constraint (|2.54|l . which is a con- 
sequence of imposing the physical requirements (|2.50() 
on the approximation l|2.46[) - l|2.48|l for non-nearest- 
neighbor interactions. Nevertheless, the form (|2.57(l sat- 
isfies Eq. (|2.50|l for arbitrary £, provided £' is still de- 
fined by Eq. (|2.55() . This strongly suggests to maintain 
Eq. H2.57[i . except that now the parameter £ is not tied to 
satisfy Eq. (|2.54|) but is instead determined from the con- 
dition y(l + ) = y(l~). Since this is not enough to guar- 
antee that y'(l + ) = y'(l~), we complement the above 
approximation by adding a linear function in the shell 
< r < 1. In summary, our low-temperature (LT) ap- 
proximation consists of 

oo 

g(r) = (r - 1)^49(1 - r) + ^ tp n {r - n)<d{r - n), (4.1) 

where the functions - 0n( r ) are given by Eq. I|2.60[) . £' 
being defined by Eq. (|2.55[) . The LT approximation l|4.1|) 
contains two free parameters (£ and A) to be determined, 
as indicated before, from the continuity conditions of y(r) 
at r = 1. To that end, note that 

y(r < 1 ) = T 4-( r - 1 ) + -|z^e-^ (4.2) 



I/(Kr<2) - L_^ e -«V + II [i + (r _i)] 
P P 
xe-^- 1 ). (4.3) 
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Now, the condition y(l + ) — y(l ) yields the following 
transcendental equation for the parameter £: 



1-x 



(4.4) 



which differs from Eq. (|2.54() in the replacement of e ^ by 
. Finally, the condition y'(l + ) = determines 
the parameter A as 



(4.5) 



This closes the construction of the LT approximation. It 
is given by Eq. (|4.1|) . supplemented by Eqs. (|2.55|) . i|2.6U[l . 
and @3J|. 
The Laplace transform of Eq. (|4.1|l is 



G{t) = A 



1-t-e- 



1 Pi(t) 

p\-Pi{ty 



(4.6) 



where the function P\(t) is given by Eq. (|2.57[1 but now 
it cannot be interpreted simply as the Laplace transform 
of the nearest-neighbor probability distribution. As al- 
ready mentioned, the LT approximation (|4.6|l fulfills the 
physical condition H2.50f) . 



r 




V(r) 




y'(r) 







i - (i - x)i'y P x 




-e^'/p i 


1 








- (1 — x)e^ ]/px 


2" 








(l-x)(l-p]/px 


2+ 


e'Ki 


- x)£' + e-t']/px 


C' 2 [l-(1 


— x)^' — e _f ]/px 



TABLE IV: Values of the cavity function and its first deriva- 
tive at r = 0, 1, and 2 in the LT approximation 14.11 . 

As an application of the explicit character of the LT 
approximation (|4.1|) , we give in Table IIVI the values of 
the cavity function and its first derivative at r = 0, 1, 
and 2. By eliminating the parameter £' in favor of y(l) 
as £' = pxy(l) one gets 



y(0) = y(l)e^« [1 - px(l ~ x)y(l)] , 



<91ny(r) 



9r 



px 2 y(l) 



r=0 



1 - - x)y(iy 



(4.7) 
(4.8) 



B. Temperature expansion 

While the approximation i|4.1[l is expected to become 
more and more reliable as the temperature decreases, it 
is not restricted a priori to low temperatures since it 
incorporates all the orders in the parameter x = 1 — x = 
e-V T ",i.e., 

Mr) = M (r) + M ] (r)3? + ^ (r)x 2 + ■■■ . (4.10) 

Let us now obtain explicitly ipn\r) and ipn\r). In the 
limit x — > 0, the solution of Eq. (|4.4I) coincides with that 

of Eq. (E3H), namely £ = £ - £ 2 e Co z H and £' = £ ~ 

£oe^°a; + - • • , where £o is defined by Eq. (|2.62(l . Therefore, 
from Eq. I2.6U[) we get 



4 0) (r) = 0, 



(0)/ 



So 



p(n-l)! 



e^V" 1 , n>l, 



(4.11) 



^)(r) = iL e -So(r-i) r n-i r^ (n + ! _ )r _ n 21 
p z n\ 

(4.12) 

From ^q 1 (r) and i/)n\r) for n > 1, and taking into ac- 
count that Eq. (|4.5|l yields A = ^e^'x 2 / p-\ — x 2 , one 

recovers the exact hard- rod result (|2.65(l . 



C. Density expansion 

In the low density regime, the solution of Eq. 14. 4|) is 

i = p{\ + xp + x*p 2 )+0{p i ) 1 (4.13) 

i 1 = xp(l + x 2 p) + 0(p 3 ). (4.14) 

Insertion into Eq. (|4.5|) yields 

A = x(l -xfp + 0{ P 2 ). (4.15) 

The density expansion of the functions ipn (f) defined by 
Eq. (|2~fJU)> is 

Mr) = (1 - x) [1 + xp{l + x-r)] + C(p 2 ), (4.16) 

^i(r) = x [l + p(x 2 + r - 2ar)] + £>(p 2 ), (4.17) 



y'(2 + )- 2/ '(2-)=p a ; 2 [y(l)] 2 . (4.9) 

Equations l|4.7|l and l|4.8|l can be considered as (approx- 
imate) extensions to finite temperatures (x ^ 1) of the 
zero separation theorems l|2.69l) and (|2.67|l . respectively. 
Equation (|4.9I) gives the discontinuity of the slope of g(r) 
at r = 2. It is interesting to note that, although Eq. 
(|4.9() has been derived from our LT approximation, it is 
not only exact in the limit x — > 1, but also in the oppo- 
site limit x — > (at px = const), as can be verified by 
comparison with Eq. Q2.36p . 



M r ) = x 2 pr + 0{p 2 ), (4.18) 

Mr) = 0(p n M, n>3. (4.19) 

Using H4.15[) - I|4.18l) in the LT approximation (|4.1I) one 
recovers the exact result y{r) = 1 + pyi(r) + 0(p 2 ), where 
?/i (r) is shown in Table ILTI Note that the inclusion of the 
term in Eq. (|4.1|) headed by the parameter A is essential 
to get this result. However, Eq. 14. 1|) does not reproduce 
the exact y2(r) for r < 2, except, of course, at x = 1. 
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D. Test for a very low temperature 
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FIG. 5: (Color online) Plot of the cavity function y(r) at 
p = 0.7 and T* = 0.15. The solid line correspond to MC 
simulations and the dashed line (which is practically indistin- 
guishable from the solid line) represents the low-temperature 
(LT) approximation 14. H . The figure also includes the exact 
(dashed-dotted line) and the PY (dotted line) cavity functions 
at zero temperature. 

As a test of the reliability of the LT approximation 
(|4.1() for very low temperatures, we present in Fig. [5| a 
comparison with MC simulations at p = 0.7 and T* — 
0.15. The agreement is excellent. In fact, the simulation 
data cannot be distinguished from the theoretical results. 
While the temperature T* = 0.15 is very low (x = 1— x ~ 
1.3 x 10~ 3 ), one can observe that y(r) in the region r < 1 
is slightly smaller than the one at T* — 0. As in Fig. 
13 Fig. [5] shows the failure of the PY approximation to 
account reasonably well for the cavity function inside the 
core. As a bonus, Fig. El illustrates that MC simulations 
of g(r) for penetrable spheres at low temperatures (say 
T* < 0.2) can be used to get (by extrapolation) the cavity 
function of hard spheres inside the core. This quantity 
is accessible from simulations of true hard spheres by an 
alternative method^ 



V. COMPARISON OF MONTE CARLO 
SIMULATIONS WITH THEORETICAL 
APPROXIMATIONS 

As Figs.0]and[S]illustrate, the HT approximation (|3.5() 
and the LT approximation (|4.1I) become very accurate 
in their expected domains of sufficiently high and low 
temperatures, respectively. However, the physically most 
interesting cases correspond to moderate temperatures 
(say 0.5 < T* < 2), where the finite penetrability of 
the particles plays a relevant role. Now some interesting 
questions are as follows: Up to what temperature the 



LT approximation remains reasonably accurate? Below 
which temperature the HT approximation ceases to be 
reliable? Do the answers to those questions depend on 
the density? To address these points we have performed 
MC simulations of the system^ for four temperatures 
(T* = 0.3, 0.8, 1.5, and 3.0) and, in each case, for three 
densities which roughly correspond to px ~ 0.3, 0.5, and 

0. 7, except for the highest temperature (T* = 3.0), in 
which case we have taken px ~ 0.5, 0.8, and 1.1. 

Figure [5] shows the simulation data and the theoreti- 
cal predictions of the radial distribution function g{r) at 
temperature T* — 0.3 and densities p = 0.3, 0.5, and 0.7. 
Although this is a rather low temperature, overlapping 
effects start to be important, especially as the density 
increases. It is then remarkable that the LT approxi- 
mation yields results that are indistinguishable from the 
simulation data. On the one hand, as one could have 
anticipated, the HT approximation does not capture the 
rich structural features that are present at this relatively 
small temperature. On the other hand, it is notewor- 
thy that, even at a temperature as low as T* = 0.3, the 
results provided by the HT approximation are not un- 
physical. 

Next we consider in Fig. the moderate temperature 
T* = 0.8 and the densities p = 0.4, 0.7, and 1. The 
LT approximation still yields excellent results, although 
it tends to slightly underestimate g(r) near contact and 
in the overlapping region as density increases. The HT 
approximation is reasonable at a qualitative level only. 

The temperature T* = 1.5 is a representative example 
of a transitional value, as illustrated in Fig.|S|for p = 0.6, 

1, and 1.4. At the lowest density (p = 0.6) the LT the- 
ory is still quite good, the HT theory being only slightly 
less accurate. At the intermediate density (p — 1), both 
theories practically coincide and underestimate g{r) near 
contact and in the overlapping region. Finally, at the 
highest density (p = 1.4), the preferable theory is the 
HT one, except for r > 1.5. Thus, although the tem- 
perature is the most important parameter to determine 
which approximation is better, the density can also play 
a relevant role. 

Figure shows g(r) at a higher temperature T* = 3 
and for the densities p = 1.7, 2.8, and 3.9. Now the 
HT approximation succeeds in reproducing the simula- 
tion data quite well for the three densities, while the 
LT approximation becomes rather poor, especially as the 
density increases. 

It is interesting to note from the simulation curves in 
Figs. EHHl that the curvature of g(r) in the shell 1 < r < 2 
changes from concave to convex as the temperature in- 
creases. Moreover, at a given temperature, the magni- 
tude of the curvature increases with the density. We 
also observe that, as temperature and density increase, 
so does the number of overlapped pairs. Nevertheless, 
the radial distribution function in the non-overlapping 
region (r > 1) becomes less structured as temperature 
increases. 

The comparison between MC data and the LT and HT 
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FIG. 6: (Color online) Plot of the radial distribution function 
g(r) at T* = 0.3 and p = 0.3 (top panel), p = 0.5 (middle 
panel), and p = 0.7 (bottom panel). The solid lines corre- 
spond to MC simulations, the dashed lines (which are practi- 
cally indistinguishable from the solid lines) represent the low- 
temperature (LT) approximation 14.11 . and the dotted lines 
represent the high-temperature (HT) approximation 13.51 . 



FIG. 7: (Color online) Plot of the radial distribution function 
g(r) at T* = 0.8 and p = 0.4 (top panel), p = 0.7 (middle 
panel), and p — 1 (bottom panel). The solid lines correspond 
to MC simulations, the dashed lines (which are practically in- 
distinguishable from the solid lines in the top and middle pan- 
els) represent the low-temperature (LT) approximation 14.11 . 
and the dotted lines represent the high-temperature (HT) ap- 
proximation 13.51 . 



theories carried out in Figs. IHHH shows that both theories 
keep being reasonably accurate well beyond their respec- 
tive expected domains. In general, the smaller the tem- 
perature and/or the density, the better the LT theory, 
while the opposite happens for the HT theory. Both the- 
ories complement each other so well that they meet and 
become practically equivalent for intermediate tempera- 



tures and densities, as illustrated by the middle panel of 
Fig- El In order to define the "basins" of both theories in 
the density-temperature plane, let us consider the locus 
of points p(T*) where the contact quantity y(l) = g(l + ) 
takes the same value in both approximations. This locus 
is plotted in Fig. ^| in the representation px vs T* . We 
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FIG. 8: (Color online) Plot of the radial distribution function 
g(r) at T* — 1.5 and p = 0.6 (top panel), p — 1 (middle 
panel), and p = 1.4 (bottom panel). The solid lines corre- 
spond to MC simulations, the dashed lines represent the low- 
temperature (LT) approximation i|4 ■ ipi ■ and the dotted lines 
represent the high-temperature (HT) approximation 13.51 . 
Note that the dashed and dashed-dotted lines in the middle 
panel are hardly distinguishable. 
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FIG. 9: (Color online) Plot of the radial distribution function 
g(r) at T* = 3 and p — 1.7 (top panel), p = 2.8 (middle 
panel), and p — 3.9 (bottom panel). The solid lines corre- 
spond to MC simulations, the dashed lines represent the low- 
temperature (LT) approximation i|4 ■ ipi ■ and the dotted lines 
represent the high-temperature (HT) approximation 13. 5H . 



have used the scaled density px rather than the true den- 
sity p because the range of values of p accessible to the 
fluid phase increases with increasing temperature roughly 
as p ~ l/xM> From Fig. 1101 and the comparison with MC 
results carried out in Figs. one can conclude that the 
LT approximation is reliable for T* < 0.8 at any den- 



sity, while the HT approximation is reliable for T* > 3 
at any density. Between these two temperatures the lo- 
cus defined by the condition ?/lt(1) = 2/ht(1) separates 
the basins of each approximation. In that transitional 
regime, the smaller (larger) the density the better the 
LT (HT) theory is. For states lying on the locus both 
approximations yield practically equivalent predictions, 
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FIG. 10: (Color online) The thick solid line represents the 
locus of points in the plane px vs T* where the LT and HT 
approximations predict the same value of the contact quan- 
tity 1/(1) = g(l + ). The LT approximation is accurate below 
the curve (shaded region), while the HT approximation is 
accurate above the curve. The circles represent the states 
considered in Figs . [^J|^| 

which slightly underestimate g(r) for distances r ^ 1.5 
(cf. Fig. |SJ. As one departs from the locus, the quality 
of the corresponding approximation (either LT or HT) 
significantly improves, as illustrated by Figs.[SH7|in the 
LT case, and by Figs. 0] and [5] in the HT case. A global 
approximation can be proposed just by adopting either 
the LT or the HT approximation, depending on which 
basin the state point lies in. More specifically, 

g(r;p,T*) = Q(a(p, T*))g LT (r; p, T*) 

+Q(-a(p,T*))g Il T(r;p,T*), (5.1) 

where 

a(p,T*)=y LT (l;p,T*)-t/HT(l;p,T*). (5.2) 

While we have been concerned in this paper with 
heuristic analytical approximations, it is worthwhile com- 
paring our MC data with the two classical integral equa- 
tion theories. For one-dimensional systems, the PY and 
HNC integral equations read 

/oo 
dsh(\ S \)y(\r- S \)f(\r-s\), (5.3) 
-OO 

/CO 
dsh(\s\) [h(\r - s\) - \ny(\r - s\)] , 
-oo 

(5.4) 

respectively, where we recall that h(r) = 
e~ tp ( r >/ kBT y(r) — 1. We have numerically solved 
the above two equations by a standard iterative method 
for the simulated states. The agreement is always quite 
good outside the core (r > 1), so that we focus now 
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FIG. 11: (Color online) Radial distribution function inside the 
core at four different states, as obtained from MC simulations 
(solid lines), the LT or HT approximations (dotted lines), 
the PY theory (dashed lines), and the HNC theory (dashed- 
dotted lines). Note that the LT results are indistinguishable 
from the MC data in the top panel, while the HNC results 
are hardly distinguishable from the MC data in the remaining 
three panels. 



on the overlapping region. Figure 1111 shows g(r) with 
r < 1 for the temperatures T* — 0.3, 0.8, 1.5, and 3, 
and, in each case, for the highest density considered in 
the simulations. Apart from the simulation data and 
the PY and HNC predictions, we have also included 
either our LT approximation (T* = 0.3 and 0.8) or our 
HT approximation (T* = 1.5 and 3), in accordance 
with Fig. Uni At the lowest temperature (T* = 0.3) 
the PY and HNC predictions seriously underestimate 
and overestimate, respectively, the simulation data, 
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while the LT approximation is excellent. However, at 
moderate and high temperatures (T* = 0.8, 1.5, and 
3) the HNC approximation becomes very good. The 
PY approximation is still worse than our LT approxi- 
mation at T* = 0.8 but is more accurate than our HT 
approximation at T* = 1.5. At T* = 3, the PY and HT 
approximations exhibit a comparable accuracy. Note 
that as temperature increases, the HT, PY, and HNC 
approximations tend to coincide and reduce to the exact 
asymptotic behavior (|2.21|) . The good performance 
of the HNC approximation at moderate temperatures 
reported here contrasts with the situation reported in 
the three-dimensional caseM^ 



VI. SUMMARY AND DISCUSSION 

It is obvious that many of the physical properties 
of real 3D systems cannot be fully represented by ID 
models. For instance, entropy-driven phase transitions, 
such as the fluid-solid transition in hard spheres or the 
demixing transition in positively nonadditive hard-sphere 
mixtures, are absent in ID. In the case of penetrable 
spheres, however, the ID model retains a number of fea- 
tures of its 2D and 3D counterparts (such as the lack 
of an exact solution, the formation of "clumps" of over- 
lapped particles*^ or the existence of a fluid-solid phase 
transition^) what justifies its study, apart from its in- 
terest at a more fundamental level. 

In this paper we have investigated the structural prop- 
erties of one-dimensional fluids of penetrable spheres. 
First we have derived the exact expression of the radial 
distribution function g(r) through second order in den- 
sity at any temperature, the results being presented in 
Table [n] This calculation already reveals some of the 
main features of the PY and HNC approximations: the 
former underestimates g(r) inside the core (r < 1), while 
the latter overestimate it. Moreover, the HNC predic- 
tions become better than the PY ones for temperatures 
larger than T* w 1. 

In the limit of asymptotically high temperatures only 
the linear chain diagrams contributing to the virial ex- 
pansion of the cavity function need to be retained. 
The series can then be resummed, the result being 
given by Eqs. (|2.21[l and l|2.27|) . The auxiliary function 
w(r), which depends on density and temperature only 
through the scaled quantity px, where x is defined by 
Eq. H2.9JI . embodies the structural properties of the high- 
temperature penetrable-rod fluid. Its series in powers of 
px converges uniformly for px < ^, so that this is the 
radius of convergence of the virial series of the equation 
of state. An explicit representation of w(r) for each shell 
n < r < n + 1 has been derived, Eq. (|2.32|) . what has al- 
lowed us to get some simple relations, Eqs. H2.35|l - H2.40[) . 

The high-temperature limit has been complemented 
by the zero-temperature limit. In that limit the system 
becomes a hard-rod fluid, whose expression for g(r) is 
known exactly. Since the particles are impenetrable in 



this limit, g(r) vanishes for r < 1. However, the cavity 
function y(r) remains finite and its determination is im- 
portant to understand the behavior of g{r) for low (but 
non-zero) temperatures. In order to derive y(r) for hard 
rods, we have first applied to finite temperatures the 
scheme providing the exact solution in the case of in- 
teractions restricted to nearest neighbors. Although the 
penetrable-rod potential at non-zero temperature allows 
for multiple interactions, the nearest-neighbor approxi- 
mation becomes more and more reliable as the temper- 
ature decreases, providing thus the exact y(r) for r < 1 
in the hard-rod limit. We are not aware of any previous 
derivation of y(r) for hard rods inside the core, although 
it might possibly exist. Quite surprisingly, the expression 
of y(r) for r < 1 is just the analytical continuation of its 
expression for 1 < r < 2, so that the hard-rod cavity 
function is fully analytical at r = 1, even though the po- 
tential is highly singular at that point. It is well known 
that the PY approximation gives the exact g(r) for hard 
rods due a fortunate mutual cancelation of the neglected 
diagrams in the region r > 1. However, it is perhaps 
less knownSi that such a cancelation does not take place 
for r < 1, what results in a wrong PY prediction of y(r) 
inside the core, already to second order in density. In 
fact, the PY and HNC cavity functions for hard rods 
exhibit an artificial second-order discontinuity at r = 1. 
The failure of the PY theory to reproduce the hard-rod 
y(r) inside the core is a precursor of its tendency to dra- 
matically underestimate the penetrability phenomenon 
at finite temperatures. 

It seems tempting to exploit the exact asymptotic be- 
havior of the correlation functions for high temperatures 
by extrapolating it to finite temperatures. Among sev- 
eral possibilities, we have taken the Pade approximant 
(|3.5|) as our high-temperature (HT) approximation. It 
is equivalent to retaining, in addition to all the linear 
chain diagrams, all the reducible diagrams factorizing 
into linear chain diagrams, although those reducible di- 
agrams have enhanced factors that partially compensate 
for the neglected irreducible diagrams. In addition, while 
x « 1/T* for high temperatures, we have kept the non- 
linear relationship (|2.9|) in the extension to finite tem- 
peratures since this is the parameter naturally appearing 
in the virial expansion. Despite its simplicity, we have 
checked by comparison with MC simulations that this 
HT approximation does quite a good job for tempera- 
tures higher than T* w 3 at any density. 

Obviously, the HT approximation is not accurate for 
low temperatures. To complement it we have constructed 
a low-temperature (LT) approximation, which is based 
on the explicit expression for g(r) obtained from the 
nearest-neighbor approximation mentioned above but in- 
troduces two changes. First, the transcendental equa- 
tion that the damping coefficient £ obeys in the nearest- 
neighbor approximation is replaced by another one com- 
ing from the continuity condition of y(r) at r = 1. Next, 
a linear function is added in the region r < 1 to en- 
force continuity of the first derivative of y(r) at r = 1. 
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More specifically, in the LT approximation g(r) is given 
by Eq. (|4.1I) . where the functions i(> n (r) are defined by Eq. 
H2.60(l and the amplitude A is given by Eq. (|4.5|) ; in those 
equations, the parameter £' is defined by Eq. (|2.55l) and 
finally the parameter £ is the solution to the transcen- 
dental equation l|4.4(l . This LT approximation reduces to 
the exact solution in the hard-rod limit. In addition, it is 
excellent for low temperatures (e.g., T* = 0.3) and quite 
good up to T* w 0.8 at any density. 

For intermediate temperatures (0.8 < T* < 3), the LT 
and HT theories compete each other, the former being 
better for densities lower than a certain threshold value 
p(T*) and the latter being better for p > p(T*). In fact, 
both theories yield almost identical predictions for states 
lying on the curve p(T*). As one departs from that curve 
by decreasing (increasing) the density and/or tempera- 
ture the quality of the LT (HT) approximation improves 
significantly. Therefore, given the simplicity and analytic 
character of both approximations, it is quite reinforcing 
that they match and complement so well that their com- 
bined use [cf. Eqs. I|5.1|) and Q5.2p ] covers satisfactorily the 
whole range of densities and temperatures. Of course, it 
would be much nicer to have a unique analytic approx- 
imation being equally accurate for both low and high 
temperatures. However, we have not been able to devise 
such a unique approach. 

To put these results in a proper perspective, we have 
also compared our MC data with numerical solutions of 
the PY and HNC integral equations. The results show 
that the PY approximation is quite poor at low temper- 
atures, slowly improving as the temperature increases. 
The HNC approximation, on the other hand, is inaccu- 
rate only at low temperatures, being excellent at moder- 
ate and high temperatures (T* > 0.8). 

Although restricted to one-dimensional systems, the 
research carried out in this paper can be used to pave 
the way to similar approaches in the more realistic case 
of three-dimensional penetrable-sphere fluids. Since the 
asymptotic high-temperature behavior is given by Eqs. 
(|2.21() and l|2.22|l for any dimensionality^^ the proposal 
of the HT approximation is straightforward. The only 
difference is that now the (three-dimensional) Fourier 
transform of the auxiliary function w(r) is given by 

(k cos k — sin k) 2 



w(k) — 96irr]x- 



(6.1) 



k 3 [k 3 — 2Arjx (k cos k — sin k)] 

where rj = (n/6)p is the packing fraction. The construc- 
tion of the LT approximation is less direct since the exact 
solution in the hard-sphere limit is not known. However, 
we can use the analytic solution of the PY equation for 



hard sphere 



,26,27.44 



as a starting point. The structure of 



the resulting proposal is 

g(r) = g(r) exp{(r - 1) [A + Br{r + 1)] 6(1 - r)}, (6.2) 
where the Laplace transform of rg(r) is 

t P(t) L -e-*(i + l + iit) 



G(t) 



12?/ 1 + P(t) ' 



P{t) 



1 + + S 2 t 2 + S 3 t 3 



In Eqs. (|6.2|l and l|6.3|l . the coefficients A, B, Li, and Si 
are functions of 77 and T* determined by imposing con- 
sistency conditions. In the limit T* — > the coefficients 
A, B, and Lq vanish and one recovers the radial distribu- 
tion function for hard spheres in the PY approximation. 
However, the cavity function y(r) in the region r < 1 
differs from the PY one even in that limit, in parallel to 
what happens in the ID case. Preliminary comparisons 
of the HT and LT approximations (|6.1|l and (|6.2|l with 
MC simulations are quite encouraging. A detailed report 
will be published elsewhere^ 
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APPENDI X: D E RIVAT ION O F EQ UATIONS 
J2.29H . J2.a0tl . AND ft'2.A'2t) 



Expanding l|2.22ll in powers of density, one gets 
with 



w n (r) 



1 

27 
1 

2^ 



dfc e lkr 



dk e ikr 



/ns(fc) 
2 sin k 



m=0 v ' J L 

(Al) 

where L is a path in the complex plane from k = — 00 
to k = +00 that goes round the singularity at k = 
from above. If r + 2m — n > 0, we can close the path 
with an upper half circle of infinite radius, so that the 
integral identically vanishes. If r + 2m — n < we close 
the path with a lower half circle, so that the contour 
encircles a pole at k = of order n. By applying the 
residue theorem, we then have 



dk- 



ik(r-\-2m — n) 



2rr 



(n-1)! 

x8(n — 2m — r) 



{-i) n {n-2m~r) n - 1 
(A.2) 



(6.3) Inserting this into Eq. (jA.ljl one gets Eq. lj2~2l?|l . 
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Let us now prove Eq. I|2.30[l . First, we note that, by 
symmetry, 



1 

2nd 



00 *k(- 2 r k)n+1 =o. 



(A.3) 



By following on the left-hand side steps similar to those 
followed to derive Eqs. IjA.ljl and (|A.2(1 . we get the iden- 
tity (for n > 2) 



[(n+l)/2] 



0- E 



(-1) 



n+m+l 



m=0 



m\(n + 1 — m)\ 



(n + l-2m)™- 1 , (A.4) 



where [a] denotes the integer part of a. Next, setting 
r = 1 in Eq. (|2.29l) and making the change m — > m + 1 , 
we have 

n x I(»+i)/2] , un+m+ i 

lu, [\ ; g 7 (n+l-2m)"- 1 . (A.5) 



m— 



m!(n + 1 — m)! 



Therefore, 

w» + i(0) 2 ro n (l) 
n + 1 n 



[(n+l)/2] n „+ TO+ l 

= (n+1) £ 



— ' mlfn. + 1 — ml! 

m=0 v ' 



x(n + 1 - 2m) r 



(A.6) 



On account of Eq. (|A.4|I . the right-hand side of Eq. I|A.6I) 
vanishes, what proves Eq. H2.30|) . 

Finally, we want to collect together all the terms con- 
tributing to the same Heaviside function in Eq. Q2.28p . 
We first define the index p = 2m — n and note that if 
p = 1 then n = 1 + 2q and m = 1 + q with q = 1,2,3,..., 
while if p > 2 then n = p + 2q and m = p + q with 
q = 0,1,2,.... Therefore, 



;(r) = E ( Xl%? q) {pXfq{l - rfqQ{1 - r) 



p=2 9 =0 

x6(p - r). 



(p + g)!g 



(A.7) 



Making use of the identity^ 

V (-^(P + ^ ^-i = 2j;(2z) , (A. 



9 =o 
7/ _ 1 



where J' p — ^(J p -i — J p +\) is the first derivative of the 
Bessel function of the first kind J p , Eq. (|A.7J| reduces to 
Eq. 
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